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Abstract 

Experiments with pedestrians could depend strongly on initial conditions. Comparisons 
of the results of such experiments require to distinguish carefully between transient state 
and steady state. In this work, a feasible algorithm - Cumulative Sum Control Chart - is 
proposed and improved to automatically detect steady states from density and speed time 
series of bottleneck experiments. The threshold of the detection parameter in the algorithm 
is calibrated using an autoregressive model. Comparing the detected steady states with 
previous manually selected ones, the modified algorithm gives more reproducible results. 
For the applications, three groups of bottleneck experiments are analysed and the steady 
states are detected. The study about pedestrian flow shows that the difference between the 
flows in all states and in steady state mainly depends on the ratio of pedestrian number 
to bottleneck width. When the ratio is higher than a critical value (approximately 115 
persons/m), the flow in all states is almost identical with the flow in steady state. Thus we 
have more possibilities to compare the flows from different experiments, especially when the 
detection of steady states is difficult. 
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1 Introduction 


In recent years, several experiments under well-controlled laboratory conditions were carried 
out to explore pedestrian characteristics in bottlenecks. Most experiments focused on the re¬ 
lationship between bottleneck width and pedestrian flow HHE]. The flow was once announced 
to grow in a stepwise manner along with lane formation [3]. However, it was shown afterwards 
by more detailed experiments to be linearly dependent on bottleneck width. The slope of the 
linear function is approximately 1.9 (m-s)“^ with bottleneck width ranging from 0.7 m to 5.0 
m isidniiis]. When bottleneck width is smaller than 0.7 m, the slope of the linear function de¬ 
creases with increasing bottleneck width [a EH HU. Other bottleneck geometrical factors were 
also studied. Shorter bottlenecks provide higher flow than longer ones laiHHii]. Wider width 
of the passage in front of bottleneck leads to higher flow mm- The flow also increases with 
increasing distance between bottleneck and holding area m- Usually bottleneck geometry was 
made up of boards higher than 2.0 m to prevent pedestrians’ bodies overlapping the boundaries. 
But Nagai et al. [1] used desks with the height of approximately 0.8 m, which was pointed out 
to actually provide wider bottleneck widths |15) . Helbing et al. HZ! and Yanagisawa et al. [I8] 
studied the influence of an obstacle in front of the bottleneck, and found that the existence 
of the obstacle leads to higher flow especially when shifted from the center. Furthermore, the 
influences of non-geometrical factors were studied. Nagai et al. [3] changed the initial density 
of the participants in the holding area, and found that the flow increases with increasing initial 
density but the rate of the increase decreases. Daamen et al. |8] changed the composition of 
the participants, and found that the experiment with mainly children has the highest flow and 
that with disabled pedestrians has lower flow. Most experiments were conducted under normal 
situations, in which the participants were asked to walk through the bottleneck with normal 
speed. Only a few experiments were conducted under competitive m, hurried [20], pushing m 
or stressful |8| situations. Unfortunately, no coincident result is made considering the influence 
of different situations. 

It is worth noting that pedestrian movement includes transient state and steady state. In 
pedestrian dynamics experiments, transient state depends strongly on initial conditions while 
steady state is a good indicator of the independency of initial conditions, especially when the 
duration of experiments is short. Thus steady state is a distinguished significance for the 
interpretation of the results, and transient state should be excluded when combining different 
experiments to get universal conclusions. However, few research about bottleneck experiments 
considered the difference between transient state and steady state. Cepolina [2T] showed the 
time series of pedestrian flow. The flow changes significantly in transient state which is at the 
beginning and at the end, but keeps constant in steady state. Yet she did not further analyse 
the flow in steady state. Seyfried et al. [6j used regression analysis to calculate the stationary 
values for the density and speed, but the flow in steady state was not considered. Rupprecht et 
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al. [inj found that the flows with all participants and with specified participants have different 
trends, but the specified participants were selected arbitrarily not considering steady state. 
Liao et al. m compared the flows in all states and in steady state, and found that they 
are both linearly dependent on bottleneck width but with different slopes. Nevertheless, the 
duration of steady state was determined according to “small fluctuations” in density and speed 
time series of the experiments, but the demarcation of the “small fluctuations” is ambiguous. 
Moreover, the specific relationship between steady state and the flow is unknown. Up to now, no 
generally accepted method was used to detect steady state in bottleneck experiments. Manually 
selected steady state can be used, but the result varies with different steady states from different 
researchers. For automatical method, Krausz et al. [22] did similar work by using a Cumulative 
Sum Control Chart (CUSUM) algorithm to detect transitions from the optical flow computations 
at the Loveparade 2010 in Duisburg. 

In this work, a feasible method to detect steady state is proposed and the steady state of 
pedestrian flow in bottleneck experiments is studied. The remainder of the paper is organized 
as follows. In section 2 a CUSUM algorithm is proposed and improved to automatically detect 
steady state, and the threshold of the detection parameter in the algorithm is calibrated using 
an autoregressive model. Section 3 applies the modified CUSUM algorithm to three groups of 
bottleneck experiments, and compares the flows in all states and in steady state. Finally, the 
conclusions are made in section 4. 


2 Detection of steady state 

2.1 Relevant variables 

In pedestrian dynamics, time series of the variables density p, speed v and flow J and their 
steady states are used to characterize the transport of pedestrian streams. The classic definition 
of flow is the number of the pedestrians passing a line per unit time: 
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J =% = 

T AU 


( 1 ) 


where At, is the time gap between the crossing of two following pedestrians, and T is the time 
for all the pedestrians N passing the measurement line. A reasonable resolution is necessary to 
identify steady state, but microscopic measurements of the flow by time gaps or mean values 
over small intervals could lead to strong fluctuations. These flow fluctuations are caused by the 
relationship between the size of the objects (here pedestrians) and the size of the facilities in 
combination with ordering phenomena like the zipper effect (see Figure 5 and Figure 7 in 0). 
To deal with this problem, we refer to the flow equation from fluid dynamics: 


J = p-vW, 


(2) 
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where W indicates the width of the facility. Here p and v is the mean density and mean speed in 
the measurement area, respectively. Following this equation, the steady states based on density 
and speed are a good indicator for the steady state of flow. 

In this paper, we aim to detect steady state from the flow in Equation (1) with relatively 
short time series. The series lengths are too short to obtain precise flow estimations, so we 
use the density and speed variables in Equation (2) to determinate if the system is in steady 
state. The aforementioned density and speed are calculated using the Voronoi method to obtain 
estimations with low fluctuations (see Equation (8) and (9) in [23]). This approach allows precise 
detection of the steady state at the scale of only few seconds. 

2.2 CUSUM algorithm 

Since the detection of steady state is a sort of transition detection, the CUSUM algorithm is 
proposed in this work. It is a sequential analysis technique which is initially used for monitoring 
transition detection |24) . The precondition of the CUSUM algorithm is a reference process with 
observations in normal situation. Let (iCi)r=i denote a sequence of observations in density or 
speed. The manually selected steady state can be regarded as the reference process 
where j > 1 and m < n (see the interval between the two green dash-dotted lines in Eigure 
1). Let Q(a) and Q(l-a) represent the upper and lower a-percentile of the distribution in the 
reference process respectively. The CUSUM algorithm continuously accumulates the 

deviations of the observations {xi)'^^^ from the a-percentiles [M] : 

sf = max{0, + Xi - Q{a)}, Sq = 0; (3) 

s~ = max{0, s~_^ Q{1 - a) - Xi}, Sq = 0. (4) 

The statistics sf and s~ show the degree of the fluctuations in (xi)^^^. If the distribution of 
(xi)^^i is different from the distribution of then the statistics increase. Therefore, the 

CUSUM algorithm is able to detect the transitions of the observations from the given reference 
process. The specific way is combining sf and s~ with a threshold of the detection parameter 
9, then the interval with both sf and s~ under 6 is the steady state. 

The preliminary CUSUM algorithm can detect the steady state, but there exist problems 
during the detection procedure. First, since the statistics sf and s~ are combined to detect 
the steady state, it is redundant to calculate them separately. Therefore, new statistics (si)f=i 
are introduced (Equation (5)). Second, the response of the statistics to the fluctuations in 
{xi)f^i is not sensitive enough. The increase rate and decrease rate of the statistics are not 
the same. To solve this problem, a step function is introduced to enhance the sensitivity of 
the response (Equation (6)). Third, the criterion in Equation (6) is not the a-percentiles of 
{xi)^j any more. The replacement is q{a), which is the upper a-percentile in a standard normal 
distribution N{0, 1). Accordingly, the standard score transformation of the observations 
is used (Equation (7)). Fourth, the maximum value of the statistics is very large that it needs a 
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Figure 1: Time series of the observations in (a) density and (b) speed. The interval 

between the two green dash-dotted lines is the manually selected steady state, which is regarded 
as the reference process where j > 1 and m < n. The interval between the two red 

dashed lines is the final steady state, which is obtained using the detected steady state minus 
the corresponding reaction time. 

long time to get back to steady state again. As an improvement, a boundary Smax is introduced 
to limit the increase of the statistics. Last but not least, since it is unknown whether at the 
beginning the statistics are in steady state or not, the value of sq should be equal to the boundary 
■Smax- From the above, the modified CUSUM algorithm is as follows: 


Si = min{max{0, Sj. 

-1 “1“ ^max}? 

^0 — ^max) 

(5) 

f 1 if 

Xj > q{a), 



■sr- 

II 


(6) 

[-1 if 

\xi\ < qia); 



Xi- n 

Xi = ; 



(7) 


where ^ is the mean of and a is the standard deviation of (xi)T^-. As shown in Figure 

2, the statistics calculated by the modified CUSUM algorithm have sensitive response 

to the fluctuations, and the trend of the statistics is clear. 

We should note that the detected transitions (the interactions in Figure 2) are not the 
real transitions in the observations (xj)^^. There is a reaction time caused by the detection 
procedure itself. When reaching steady state and leaving steady state, the reaction time is 
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(a) density (b) speed 

Figure 2: Detection of steady state by the modified CUSUM algorithm. The statistics in (a) 
density and (b) speed are calculated based on the observations in Figure 1 (a) and (b), 

respectively. Si shows the degree of the fluctnations in the observations The interval 

with Si under 6 is the detected steady state. 


calculated as follows, respectively: 


^reaching — 


-6 


f 




eaving — 


e 

7 ’ 


( 8 ) 


where / is the frame number per second in the observations. The hnal steady state is using the 
detected transitions minus the corresponding reaction times (see the interval between the two 
red dashed lines in Figure 1). 


2.3 Threshold of detection parameter 

The accurate detection of steady state depends on the detection parameter’s threshold 9 to be 
calibrated. For a Neyman-Pearson statistical test, 0 ( 7 ) is the upper 7 -percentile of the statistics 
for the reference process In this way, we are able to control the probability of 

detecting a false transition, which is the risk level I- 7 . 

Different methods to calibrate 9 are shown in Figure 3. The simplest way is the direct 
calibration with the reference process {j > 1 and m < n). This method requires that 

the number of the selected observations in the reference process should be large enough for 
calibration. If not, a bootstrap method is used to expand the reference process [22]. The 
reqnirement is that the observations in {xi)^j should be independent of each other, or that the 
number of the observations in (xi)^^ should be large enough (moving block bootstrap). 
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Figure 3: Different methods to calibrate the threshold of the detection parameter 0. 


In this work, the observations do not meet any requirements of the two empirical methods 
in Figure 3. Therefore, an autoregressive process is proposed to model the standard 

score of the reference process 


Vi = c-Vi 


_i + 


2/0 = 0 , 


(9) 


where c is the first autocorrelation of and (ei)^i are independent normal random vari¬ 

ables. The autoregressive process {yi)f^i is stationary with normal distribution, and cor(yj, Ui+i) - 
c for all i. Using the autoregressive model, the threshold 0 ( 7 ) can be obtained by simulations 
or analytically by using mathematical tools of stochastic processes. If by simulations, the simu¬ 
lation time T should be sufficiently long enough to exclude fluctuations and precisely estimate 
the distribution of Moreover, several runs must be done to control the precision of the 

estimation. When c, 7 ~ 1, the computational effort to get good estimations is important for 
standard computers. For instance if c = 7 = 0.99, the process has to be simulated T = leS steps 
to obtain precise results. The threshold 0 ( 7 ) with the autoregressive process can be faster ob¬ 
tained analytically. The couple (y*, Si) is a Markov chain with stationary distribution |u(x, s) dx 
(x G M, 0 < s < Smax) such that (global balance equation): 


fi{x,s) = 


j I y{y,s-l)g{y,x)dy if \x\ > q{a), 


\j y{y,s + l)g{y,x)dy if \x\<q{a), 
for all 0 < s < Smax) and at the borders 

/^(a;,Smax) = y'(/i(2/,'Smax - 1)+h(2/,'Smax))5'(2/,a:)dy and g{x,0) = 0, 

if |x| > q{a), and 

g{x,0) = j{g{y,l) + iJ,{y,0)) g{y,x)dy and y(x,Smax) = 0, 


( 10 ) 


( 11 ) 


( 12 ) 
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if |x| < q{a). Here g{y, x) = ^ 75 ;^ exp • 

Unfortunately, Equations (10-12) do not have an explicit solution for g{- ,•). Yet they can 
be approximated by using the rectangular numerical scheme Xi = —x' + Sxi with Sk = 2 x'/ K, 
ai,s = fJ'ixi, s) and bi^k = Sk g{xk, Xi). The numerical approximation yields in the linear equation 


Mx = 0, 


(13) 


where 


X = "^(ao.o, ai,o, • • • aK,o, ao,i, ■ ■ ■ aK,sm^^) 


(14) 


and M = B — Id with 
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0 
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bK,0 • • 

■ bK,K 





(15) 


ia = argmaxj{xi < —q{a)} and % = argminj{xj > q{a)}. The stationary distribution of the 
couple {ni, Si) is approximated by solving a linear system with {K + l)(smax +1) equations. The 
complexity of the resolution is in 0{sraa.xK^) by using the Thomas-Algorithm (see in Appendix), 
which makes in general the numerical approximation of the analytical solution faster than the 
simulations. It is difficult to estimate the error of the numerical approximation. Tests show 
that A' ~ 50 gives good results. The error can be reduced by using trapezoidal or polynomial 
schemes instead of rectangular ones. 


2.4 Robustness 

For a series of observations (xj)(U]^, the steady state detected by the modified CUSUM algorithm 
depends on the combination of the statistics (si)”=i and the threshold of the detection parameter 
6. Scrutinizing Equations (5-7) and (9), six key parameters determining the values of (si)iU^ 
and 9 are summarized in Table 1. a, 7 and Smax are independent of the observations so they 
should be calibrated independently. The value of a should be close to 1 to diminish outliers in 
(xi)”ij as much as possible. The value of 7 should be close to 1 to diminish false detections 
as much as possible. The boundary Smax of the statistics ('Si)(Li should be approximately two 
times of 6. In this work, the values of the three calibrated parameters are fixed, and a = 0.99, 7 
= 0.99, Smax = 100. On the contrary, //, a and c vary with different observations because they 
are measured from the original reference process {xi)^y g, and a affect the statistics (si)(Li 
(Equation (7)), and c affects the threshold of the detection parameter 6 (Equation (9)). In the 
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Table 1: Key parameters in the modified CUSUM algorithm. 



a 

critical probability for the reference distribution 

Calibrated parameters 

7 

confident level 


^max 

boundary of the statistics (si)f=i 


u 

mean of the observations in {xi)^j 

Measured parameters 

a 

standard deviation of the observations in {xi)^j 


c 

autocorrelation of the observations in (xi)T^- 


situation that and 0 both change with different reference processes, the detected steady 

states might be different. 

To investigate the robustness of the modified CUSUM algorithm, three reference processes 
are manually selected from the same observations and they do not overlap with each other 
(see Figure 4). Keeping all the conditions the same, the detection procedure is repeated for 
each reference process. The results show that the three detected steady states in Figure 4 are 
almost the same. Moreover, they coincide with the steady state in Figure 1 (a). Therefore, for 
a given series of observations, the detected steady state is fixed with any reasonable reference 
process from the observations. By using the modified CUSUM algorithm, researchers are able 
to reproduce the same steady state with different manually selected reference processes. 

3 Pedestrian flow in bottleneck experiments 

3.1 Experimental setup 

Three groups of bottleneck experiments are analysed in this work. The principle setup is shown 
in Figure 5. Bottleneck width and bottleneck length is h and I, respectively. The width of the 
passage in front of the bottleneck is w. The distance between the bottleneck and the holding 
area is d. Before each run, N participants are arranged in the holding area with initial density 
Pini. The corresponding parameters in each experiment are listed in Table 2. The experiment 
EG was conducted in 2006 in Diisseldorf, and the experiments AO and UO were conducted 
in 2009 in Diisseldorf. More information of the experiment EG, AO and UO is given in |10j . 
nais] and [23], respectively. Note that the shape of the holding area in the experiment AO 
was not rectangular but semi-circular with radius 8.6 m directly in front of the bottleneck (see 
Figure 1 in |15ji. Also note that d in the experiment UO was 12.0 m, which included a free 
region of 4.0 m long and a passage of 8.0 m long (see Figure 1 (a) and (b) in |53|). 

The modified CUSUM algorithm is applied to detect the steady states of the three groups of 
bottleneck experiments. For the experiments EG and AO the measurement area of the detected 
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Figure 4: Time series of the observations in density. The red slashed, green latticed 

and magenta back-slashed areas represent the manually selected reference processes. The cor¬ 
responding detected steady state is the interval between the two red dashed lines, the two green 
dash-dotted lines and the two magenta dotted lines, respectively. 

variables, density and speed, is at the beginning of the bottleneck with size b (see the region 
slashed in green in Figure 5). For the experiment UO the measurement area is in front of the 
bottleneck with size b (see the region back-shaded in blue in Figure 5). The measurement line 
of pedestrian flow is always at the beginning of the bottleneck. 

3.2 Steady state and flow 

The relationship between flow and bottleneck width is shown in Figure 6. All the flows are 
calculated according to Equation (2). The error bars represent the standard errors of the flows 
in steady states. The calculation of the standard error is by dividing the standard deviation 
by the square root of the number of the measurements that make up the mean. Since the 
observation (flow) in each frame is dependent with the observations nearby, they cannot be 
used as the measurements directly. Thus the observations in steady state are divided into 
groups that the number of the groups P is the number of the measurements to calculate the 
standard error: 

P = -, (16) 

P 
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Figure 5: Experimental setup, b is bottleneck width. I is bottleneck length, w is the width of 
the passage in front of the bottleneck, d is the distance between the bottleneck and the holding 
area. Before each run, N participants are arranged in the holding area with initial density /?ini. 
The measurement area of density and speed for the experiments EG and TO is slashed in green, 
and that for the experiment UO is back-slashed in blue. The size of the measurement areas is 
b*l. 

where F is the number of the observations in steady state, and p is the number of the obser¬ 
vations in each group. Here p should be large enough to make the groups independent of each 
other, but also small enough to ensure the sufficient number of the measurements to diminish 
fluctuations when making up the mean. In this work, p = f and if P < 20 the bias of an 
estimator is used. As shown in Figure 6, the flow is linearly dependent on bottleneck width but 
the slope of the linear function varies in different experiments. This reconfirms the result that 
the flow per unit width is constant as bottleneck width changes mm- 

For each run in each experiment, the flows in all states and in steady state are compared 
and the difference between them is analysed. As shown in Figure 6, the difference changes 
significantly in the last few runs of the experiments EG and AO but slightly in the experiment 
UO. To scrutinize what exactly impacts the difference, the parameters in Table 2 are studied. 
Bottleneck length I, passage width w, the distance between bottleneck and holding area d, 
pedestrian number N and initial density pi^i are same in the experiments EG and AO. Since 
the difference changes significantly in part runs of these two experiments, the above parameters 
have no direct influence on the difference. Bottleneck width b is positively correlative with the 
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Table 2: Parameters in three different bottleneck experiments EG, AO and UO. 6 is bottleneck 
width. I is bottleneck length, w is the width of the passage in front of the bottleneck, d is the 
distance between the bottleneck and the holding area. N is the number of participants, pini is 
the initial density of the participants in the holding area. 


Experiment 

EG 

AO 

UO 

b [m] 

0.9, 1.0, 1.1, 1.2, 1.4 

1.6, 1.8, 2.0, 2.2, 2.5 

2.4, 3.0, 3.6, 4.4, 5.0 

0.7, 0.95, 1.2 

1.0, 1.3, 1.6 

0 .8, 1.2, 1.6, 2.0 

1 [m] 

4.0 

1.0 

0.1 

w [m] 

7.0 

18.0 

1.8, 2.4, 3.0 

d [m] 

4.0 

0 

12.0 

N [persons] 

180 

350 

150, 250, 400 

Pini [m“2] 

2.6 

3.0 

3.0 

N/b [persons/m] 

200, 180, 164, 150, 129 

113, 100, 90, 82, 72 

146, 117, 97, 80, 70 

214, 158, 125 

250, 192, 156 

500, 333, 250, 200 


difference in the experiments EG and AO, but the relationship is ambiguous in the experiment 
UO. In this situation, we cannot judge that bottleneck width b is directly proportional to the 
difference. Afterwards, different combinations of the parameters in Table 2 are studied. The 
ratio of pedestrian number to bottleneck width N/b is found to have a direct influence on the 
aforementioned difference, which is represented by Z in Figure 7. Z is calculated by using the 
flow in all states minus the flow in steady state and taking the absolute value. As shown in 
Figure 7, Z decreases with the increase of the ratio N/b. The critical value of the ratio N/b is 
approximately 115 persons/m, which determines if the difference between the flows in all states 
and in steady state is acceptable. Considering the practical significance of N/b, it describes 
whether the number of participants is large enough for the corresponding bottleneck width to 
reach a steady state. When the value of N/b is higher than the critical value, the steady state 
is quickly reached and the duration of steady state is long. Conversely, it takes a long time to 
reach steady state and the duration of steady state is short when the value of N/b is lower than 
the critical value. In the case of extremely small values of N/b, the steady state even cannot 
be reached. From the above, the ratio of pedestrian number to bottleneck width N/b is the 
main influence factor of the difference between the flows in all states and in steady state. Other 
factors, such as the composition of participants and the motivation of the experiments, also 
impact the difference aforementioned, but they are beyond the consideration of this work. 
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Figure 6: Relationship between flow J and bottleneck width b. For each experiment, the flows 
in all states and in steady state are represented by same symbol with different colours. The 
error bars represent the standard errors of the flows in steady states. The black dashed line 
corresponds to b = 0.4 m, which is the minimum width for pedestrians to pass. 

4 Conclusions 

In this work, a feasible CUSUM algorithm is proposed to automatically detect steady state 
from density and speed time series of bottleneck experiments. To improve the algorithm, a step 
function is introduced to calculate the statistics for enhancing the sensitivity of the response 
to the fluctuations. In addition, a boundary is added to limit the increase of the statistics. 
The threshold of the detection parameter in the algorithm is calibrated using an autoregressive 
model. Comparing the detected steady states with previous manually selected ones, the modified 
CUSUM algorithm gives more reproducible results. 

For the applications, three groups of bottleneck experiments are analysed by the modified 
CUSUM algorithm. The steady states of density and speed in each run are detected separately. 
Then the interval which is included both in the steady states of density and speed is regarded 
as the steady state of the flow. The flows in all states and in steady state are measured and 
compared. The results reconfirm that the flow per unit width is a constant as bottleneck width 
changes. Furthermore, the difference between the flows in all states and in steady state mainly 
depends on the ratio of pedestrian number to bottleneck width. The critical value of the ratio 
is approximately 115 persons/m. When the value of the ratio is higher than the critical value. 
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Figure 7: Relationship between the difference Z and the ratio N/b. Z is calculated by using 
the flow in all states minus the flow in steady state and taking the absolute value. The black 
dashed line corresponds to N/b = 115 persons/m, which is the critical value of the ratio N/b 
to determine if the difference Z is acceptable. 

the steady state is quickly reached and lasts longer. In this situation, the flow in all states is 
almost identical with the flow in steady state. Considering the ratio, we have more possibilities 
to compare the flows from different experiments, especially when the detection of steady state 
is difficult. 

In future studies, the flow in steady state should be used when combing the flows from 
different bottleneck experiments. For the experiments in which it is difficult or impossible to 
detect steady state, the ratio of pedestrian number to bottleneck width should be calculated to 
estimate the difference between the flows in all states and in steady state. 
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Appendix 

The matrix M in Equation (13) has a size of {K + l)(smax + 1) x {K + l)(smax + !)■ Solving M 
directly with Gaussian elimination takes 0{s^^^K^) operations, which is not feasible for great 

^max* 

Since M is of block-tridiagonal structure, a Block Thomas-Algorithm [25] is used to solve 
the system efficiently in O(sniaxAr^). Essentially it is a block matrix variant of the Thomas- 
Algorithm, which is a simplified version of the Gaussian elimination on tridiagonal matrices. 
The outline of the algorithm is as follows: 


1 

o 


Xl 


ri 

Li D2 U2 


X2 


r2 

L2 D3 



= 


Us 

*max 


Xsmax 


"^Smax 

0 Lc Dc -1-1 

Amax Amaxu^-L 


^Smax + l 


^Smax-l-1 


where Di, Ui,Li are matrices of size (K -|- 1) x (K + 1) and Xi, Vi are vectors of size K + 1. 
Since the zero blocks on the off diagonals are not necessary, using a sparse storage scheme 
reduces the memory amount for the matrix from ((AT -|- l)(smax + 1))2 to (3 '^max + 1){K + If. 
To obtain a solution, modify the entries in a forward sweep: 

Ul = 

r* = 

for i = 1, and 

U* = (A - LiUUr^Ui, 

Tl = {D,-LiUU)-\n-UTU), 


for % — 2, 3,, Sjiiax T !• 

Then in a backward sweep: 


^Smax-l-l ^Smax + 1’ 


for i = Smax + 1, and 


Xi = r* - U*Xi+i, 


for i = 


’max? '^max 


- 1 , 


, 2 , 1 . 
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